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I review recent developments in determining the QCD phase diagram by means of lattice simu- 
lations. Since the invention of methods to side-step the sign problem a few years ago, a number 
of additional variants have been proposed, and progress has been made towards understanding 
some of the systematics involved. All available techniques agree on the transition temperature 
as a function of density in the regime li q /T < 1, There are by now four calculations with signals 
for a critical point, two of them at similar parameter values and with consistent results. However, 
it also emerges that the location of the critical point is exceedingly quark mass sensitive. At the 
same time sizeable finite volume, cut-off and step size effects have been uncovered, demanding 
additional investigations with exact algorithms on larger and finer lattices before quantitative con- 
clusions can be drawn. Depending on the sign of these corrections, there is ample room for the 
eventual phase diagram to look as expected or also quite different, with no critical point at all. 
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Figure 1: Qualitative QCD phase diagram for Nf = 2 according to general expectations. For Nf = 3 and 
m < m c , the diagram looks as on the left with the transition being first order all the way, while for m > m c it 
looks as on the right. 



1. Introduction 

1.1 Qualitative expectations 

In the physics communities dealing with QCD at finite temperature and density, there appears 
to be little doubt that the (T,fi) phase diagram qualitatively looks as in Fig. |l| (right). Given that 
until recently non-perturbative calculations were impossible for /I ^ 0, and even on the temperature 
axis simulations with dynamical fermions have only become feasible in the last few years, it seems 
worthwhile to recall the qualitative arguments that lead to this picture. Such a diagram represents 
one set of quark masses. As theorists, we also view the quark masses as parameters and wish to 
understand the phase diagram of the entire parameter space {m u j,m s ,T,iJ.}, which should aid us 
in unveiling the physical situation as well. 

Starting point of the argument §\§ is the Nf = 2 theory with degenerate quark masses. At zero 
density and in the chiral limit, pt , m = 0, the chiral condensate represents a true order parameter 
distinguishing between separate phases, and the symmetry breaking pattern is SU (2)y x SU (2) a — > 
SU (2)y. A local order parameter vanishing everywhere in one phase and being non-zero in another 
corresponds to a non-analytical function of the parameters of the theory, thus requiring a true phase 
transition and excluding an analytical crossover. If the corresponding phase transition is second 
order, then chiral symmetry implies that it should be in the universality class of 3d 0(4) spin 
models, a scenario which has been very popular among theorists. Note, however, that a first order 
transition is a logical possibility as well. 

For low temperatures and large densities, a number of model calculations (see e.g. [§]) ap- 
pear to agree on a first order transition between nuclear matter and quark matter in a colour- 
superconducting state, which is expected for asymptotically high density. Fig. [T] represents only 
the simplest picture, other variants have one or several phases between the hadronic and the super- 
conducting phase, see [Q, ||] for more details. The most natural scenario then has the first order 
line at finite density joining up with the second order line coming from the temperature axis in a 
tri-critical point. On the other hand, if the quarks have finite masses, chiral symmetry is explicitly 
broken and there is no true order parameter. The chiral condensate still experiences a rapid change 
of value at the pseudo-critical temperature, but now it is an analytic crossover. In this case the first 
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Figure 2: Left: Pseudo-critical coupling and temperature defined by the peak of a susceptibility. Middle: 
Schematic phase diagram for Nf = 3. Right: Projection onto the critical surface. 

order line at finite density has to terminate in a critical endpoint. 

For three degenerate flavours Nf = 3, pL = 0, the chiral limit exhibits a first order phase tran- 
sition. First order transitions are stable under small variations of the parameters, and thus the first 
order regime extends to small masses m < m c , for which it can actually be measured on the lattice. 
With increasing quark mass the transition weakens until it ends in a critical point at m c , and for 
m > m c a smooth crossover is observed. For the (r,ju) phase diagram this implies a first oder line 
connecting the transitions on the axes for m <m c , while for m > m c the first order line emanating 
from the /i-axis again has to terminate in a critical endpoint. 

The "standard scenario" for the physical case with Nf = 2+1 quarks is as in Fig. [I] (right), with 
the critical endpoint moving to larger pL c with increasing quark masses, as may be inferred from a 
continuity argument. For Nf = 3 with m <m c the phase diagram has a first order line connecting 
both axes. Upon sending m s — > °° this picture should continuously evolve into the Nf = 2 diagram 
with a critical endpoint, thus implying d\l c jdm s > 0. 

The full phase diagram of the Nf = 2,3 theories is in the 3d space {m, T,[j,}, as in Fig. ^| In 
order to map it out by simulations the first step is to identify the critical surface 7o(/i,m) separating 
the high and low temperature regions. Since simulations are always on finite volumes, this surface 
is only pseudo-critical and represents a smooth crossover. It can be defined by, e.g., peaks in 
susceptibilites, cf. Fig. ^ (left). This step is typically rather straightforward. The much more 
difficult task is to perform a finite size scaling analysis to identify the order of the transition in the 
infinite volume limit for the different regions of parameter space. For Nf = 3, such an analysis 
yields a critical line separating a first order region from a crossover region on the surface To(/J,,m), 
Fig. U (middle). It is convenient to eliminate the temperature axis from this diagram by projecting 
onto the pseudo-critical surface, i.e. temperature is always implied to be To(n,m), Fig. ^| (right). 

This form of a phase diagram is particularly suitable to display the three flavour theory with 
non-degenerate masses, Nf = 2 + 1, including the special cases Nf = 2,3, as in Fig. ||. Note that 
because of the difficulties of simulating dynamical fermions, even for pL = (left) we know very 
little about the critical lines separating the first order from the crossover regions. Up to now the 
only published point that has been calculated to some accuracy with standard staggered fermions 
is the critical point m u ^ = m s = m c on the Nf = 3 diagonal [ffl, BL [6|], which was numerically 
identified to belong to the universality class of the 3d Ising model [Q]. While the statement about 
the universality class concerns infrared physics and thus is stable against cut-off effects, the location 
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Figure 3: Left: Schematic phase diagram for Nf = 2+1 at jx = 0. Temperature is implied to be (pseudo-) 
critical, To(m u c j,m s ), everywhere. Right: The same with finite quark density. 

of the critical point in the bare mass diagram is very sensitive to renormalisation effects. To date 
only N t = 4 calculations (a ~ 0.3 fm) have been performed, but simulations with improved actions 
give values for m c which are about ~ 1/4 of the standard action result [Q]. The critical line for 
non-degenerate quark masses is being calculated presently, cf. |{7]] and Section All available 
results are consistent with the physical point lying on the crossover side of the boundary. This has 
also been found in a recent simulation with standard staggered quarks with a pion to rho mass ratio 
tuned to its physical value [|8]J. 

When a finite quark number density is switched on, a /x-axis for the chemical potential has to 
be added to the diagram, and the critical line separating the first order region from the crossover 
region turns into a critical surface, as indicated in Fig. |3] (right). The standard scenario with m c (ji) 
being an increasing function of the chemical potential then implies that this surface bends towards 
larger quark masses. Consequently, tuning the quark masses to the physical point and switching 
on a chemical potential, the intersection with the critical surface marks the critical value juv of the 
end point, beyond which there is a first order transition. Thus, a determination of the QCD phase 
diagram in the full parameter space {m u> d,m s , T,n} entails mapping out these critical surfaces and 
understanding how they are joining up in the different limit theories. 

1.2 Lattice QCD at finite temperature and density 

Standard Monte Carlo simulations at finite density are made impossible by the so-called sign 
problem of the lattice grand canonical partition function, 

Z = J DUDxfrD\ife~ s ^~ s f [u ^ = J DU [detM(ii)] f e~ s * [u] , S f = Y,W¥- (1-1) 

For /i = 0, the relation 75M75 = M 1 " guarantees positivity of the fermion determinant, detM > 0, in 
every gauge background. For the gauge group SU(3), the fermion determinant becomes complex as 
soon as a non-zero quark chemical potential /J. = /Ib/3 is switched on. Thus it cannot be interpreted 
as a probability distribution, which rules out standard importance sampling. 
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The sign problem of QCD is still unsolved to date. Successful solutions of fermion sign 
problems in a number of spin models by means of cluster algorithms ^ unfortunately do not seem 
to generalize to QCD. However, significant progress has been made since 2001 with a number of 
different approaches that circumvent the sign problem, rather than solving it. These approaches 
can be put into three categories: 

• Two-parameter reweighting 

• Taylor expansion in \ijT 

• Simulations at imaginary \i, either analytically continued to real \i or Fourier transformed to 
the canonical ensemble 

All these methods have limitations and presently work reliably only for small enough fi/T<l. 
However, the systematics is different between them, thus allowing for meaningful cross checks 



(for a comparison of early results, see [ 10]). Another alternative is to study related theories without 
the sign problem, such as SU(2) QCD and QCD at finite isospin, the latter being close enough to 
the case of interest for meaningful comparisons. The fact that all agree in the determination of the 
pseudo-critical temperature To(mi,ix) is one reason for the recent enthusiasm in this field, and gives 
reason to hope that the order of the transition may be settled in the near future as well. 

2. Massless Nf = 2 at zero density: 0(4) or first order? 

Before delving into the discussion of finite density calculations, let us turn to the /I = be- 
haviour of the theory which played an important role in the derivation of the qualitative phase 



diagram in Section 1.1. It is a longstanding question whether the phase transition in the chiral limit 



of the two-flavour theory is indeed second order with 0(4) universality or first order. On the lattice, 



0(4) will effectively look like 0(2) as long as there are discretisation effects [11]. A lot of work 



has been done over the years, but no definite conclusion has been reached. Among the more recent 



work, Wilson fermions appear to see 0(4) scaling Q12[], while staggered actions are inconsistent 
with both 0(4) and 0(2) [|l3|]. (The staggered strong coupling limit, however, does display 0(2) 
scaling [0]). 

A new attempt to tackle this question by means of a finite size scaling analysis with unprece- 
dented lattice sizes was made in [15]. The work simulates L 3 x 4 lattices with L = 16 — 32, using 



the standard staggered action and the hybrid Monte Carlo R-algorithm [|16|]. Several quark masses 
are studied, the smallest being m/T< 0.055. In a critical region quantities like, e.g., the specific 
heat or the chiral susceptibility scale universally as 



C v - Co ~ L a/V f c ( %l}l v ,amV h \ , T = 1 — T/T c 

X l?> v f x (%l}l v ,amiy^ . (2.1) 
Here the non-singular part of the specific heat Co has been subtracted. The values for the exponent 



yh are known with some precision and nearly the same for 0(4) and 0(2). The authors of [15] 
thus fix yn to this value, and then choose L and m for a series of simulations such as to keep 
(amL yh ) constant. This reduces the two-parameter scaling problem to depend on one remaining 
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Figure 4: Finite volume scaling behaviour with of specific heat and chiral susceptibility. For 0(4) viz. 0(2) 
behaviour, the data should fall on a horizontal line [n5D . 
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Figure 5: Testing the mass scaling of the chiral susceptibility for 0(4) and first order behaviour [JT3]. 



variable only, which can be more easily scanned. The infinite volume limit in this procedure thus 
corresponds to the chiral limit and allows to check whether the data are consistent with the predicted 
scaling behaviour. 

Fig. || shows simulation results from [15]. Scaling as in Eq. ( [01 ) would imply the data points to 
fall on a horizontal line, which is clearly not the case. Alternatively, one may keep the other scaling 
variable rL 1//y fixed and vary the quark mass. Furthermore, in place of 0(4) or 0(2) exponents, 
also consistency with first order exponents can be tried. Some results of this attempt are shown 
in Fig. |5[ The fit to first order scaling is slightly better, but not very convincing either. Moreover, 
D'Elia et al. looked in detail at plaquette distributions in the transition region as well as Monte 
Carlo histories. No signs of a metastability region commensurate with a first order transition could 
be observed. 



Many variations of the analysis are performed in [|15J, with similar outcome. Hence, even after 
formidable computational effort the question of 0(4) vs. first order scaling remains open for the 
moment. There are several possible explanations. The scaling region away from the chiral limit 
could be exceedingly small, or discretisation effects could play a large role. This is also suggested 
by the fact that Wilson and staggered fermions appear to scale differently. The question of cut-off 
effects will be addressed by the authors of [ 15], who announced an investigation at N t = 6 as well. 
Another possibility is that exceedingly large volumes might be required to distinguish weakly first 
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order and crossover. An example for such behaviour is two-colour QCD in the strong coupling 
limit. Numerical results for this theory require lattice sizes L > 128 before the correct scaling is 



observed [|17[]. An important observation is that, for both scenarios in Fig. g, it is the lowest mass 
data points which spoil the fits. This indicates possible systematic errors of the Monte Carlo for 
very low masses. We shall indeed see in Section ^3] that in this regime the R-algorithm has strong 



step size effects for steps of half the quark mass, as chosen for the lowest mass point in [|15|]. These 
effects change the apparent order of the phase transition. Thus, for any future investigation an exact 
algorithm is necessary. Meanwhile, we should keep an open mind to the possibility of a first order 
transition in the chiral limit. In this case the phase diagram Fig. [j] would be as for Nf = 3 with a 
very small critical quark mass m c . There would be a first order line all the way for m < m c , or a 
first order line with an endpoint, Fig. [j] (right), for m > m c . 

3. Finite density phase diagram from two parameter reweighting 

Significant progress enabling finite density simulations was made a few years ago, by a gener- 



alisation of the Glasgow method [ 1 8 ] to reweighting in two parameters [19]. The partition function 
is rewritten identically as 

rw _l e-^)det(MGu)) \ 



e- s ,(A.)det(M(M = 0))/ M=OA) 



where the ensemble average is now generated at /I = and a lattice gauge coupling jS , while a 
reweighting factor takes us to the values pL , j8 of interest. The original Glasgow method reweighted 
in /I only and was suffering from the overlap problem: while the reweighting formula is exact, 
its Monte Carlo evaluation is not. The integral gets approximated by a finite number of the most 
dominant configurations, which are different for the reweighted and the original ensemble, and this 
difference grows with }X. When calculating critical behavior at some pL, one-parameter reweight- 
ing uses a non-critical ensemble at ju = 0, thus missing important dynamics. By contrast, two- 
parameter reweighting proceeds along the pseudo-critical line of the phase change, thus always 
working with an ensemble that probes both phases. This approach produced the first finite den- 
sity phase diagram from the lattice, obtained for light quarks corresponding to m n ~ 300 MeV 



A Lee- Yang zero analysis [ |21| ] was employed in order to find the change from crossover 
behaviour at \i = to a first order transition for \i > pL c . A later simulation at physical quark 
masses puts the critical point at fig ~ 360 MeV [g], Fig. ^ (left). In this work L 3 x 4 lattices 
with L = 6 — 12 were used, working with the standard staggered fermion action and using the R- 
algorithm. Quark masses were tuned to m u j/To 0.037, ra. s /7o « 1, corresponding to the mass 
ratios m n /m p « 0.19, m n jm^ f» 0.27, which are close to their physical values. 

A difficulty in this approach is that the determinant needs to be evaluated exactly. Because 
of the sign problem the reweighting factor is exponentially suppressed with volume and chemical 
potential, thus limiting the applicability to moderate values of those parameters. This point will be 
discussed in more detail later. Moreover, since the statistical fluctuations are those of the simulated 
ensemble instead of the physical one, it remains difficult to obtain reliable error estimates. For a 
proposed procedure see [|22|]. 
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Figure 6: Left: The phase diagram for physical quark masses as predicted by the two parameter reweighting 
method [ph. Right: High density results from the density of states method, indicating a triple point | 



In present work in progress two-parameter reweighting is combined with the density of states 
method [23], in order to extend the applicability of reweighting to larger values of ji/T and thus 
to lower jl. First interesting results, with indications for a possible triple point, are shown in Fig. || 
and presented in more detail in these proceedings [24]. 



4. Finite density by Taylor expansion 

Another method to gain information about non-zero pL is to compute the coefficients of a Taylor 
series expansion of observables in powers of /J./T. Early attempts have looked at susceptibilities 
and the response of screening masses to chemical potential [^, 26, 27, 28]. More recently it has 
also been used to gain information on the phase transition and its nature itself [29]-[32]. This idea 
exploits the fact that on finite volumes there are no non-analytic transitions, and hence the partition 
function Z(m > 0,/J., T) is an analytic function of the parameters of the theory. For small enough 
H/T one may then hope to get away with only a few terms, whose coefficients are calculated at 
/i = 0. Moreover, CP symmetry of the QCD action translates into a reflection symmetry of the 
partition function, Z(jx) = Z(— ju), such that real physical observables have series expansions in 
(/i/r) 2 . Thus, in particular the pressure density can be expressed as an even power series, 

p (r, M) = -£ = (^iogz(r,M), ^ = t c 2n(T)^) 2n . (4.i) 

Since only even terms appear, the coefficients are equivalent to generalised quark number suscep- 
tibilities at ji = 0, and hence measureable with standard simulation techniques. For high enough 
temperatures T > Tq, the scale of the finite temperature problem is set by the Matsubara mode 
~ 7lT, and one would expect coefficients of order one for an expansion in the 'natural' parameter 
jj. I (%T) [H]. We shall see later that this is borne out by simulation results. 
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Since all the jj, -dependence in the partition function is sitting in the fermion determinant, it is 
derivatives of the quark matrix that need to be computed, 

dlndetM / ,dM\ dtrAr 1 / , dM A 

tr I AT" -5— ] , = -tt > etc., (4.2) 



d/jL \ djJ, J ' d/x \ dpi 

which can be iterated for higher orders. These expressions become increasingly complex and 



methods to automatize their generation have been devised Q32p . Note that one now is dealing 
with traces of composite local operators, which greatly facilitates the numerical evaluation in a 
simulation compared to a computation of the full determinant. The numerical estimate of these 
expressions proceeds by the random noise method, with typically O(10) — 0(100) Gaussian noise 
vectors. 

If one is interested in phase transitions, finite volume scaling towards the thermodynamic limit 
has to be considered. True phase transitions will emerge as non-analyticities in the pressure, which 
is not the case for analytic crossover behaviour. Given that in the two flavour theory with finite 
masses and for /j. = the deconfinement transition is an analytic crossover, one may expand about 
/i = and then look for the emergence of a finite radius of convergence as the volume increases. 
The radius of convergence of a power series gives the distance between the expansion point and the 
nearest singularity, and may be extracted from the high order behaviour of the series. Two possible 
definitions are 



p,r = lim p n ,r n with p„ 



Cln 



1 /In 



Cln 



Cln+1 



1/2 

(4.3) 



General theorems ensure that if the limit exists and asymptotically all coefficients of the series are 
positive, then there is a singularity on the real axis. More details as well as previous applications 
to strong coupling expansions in various spin models can be found in [03p. In the series for the 
pressure such a singularity would correspond to the critical point in the (ju, T) -plane. 

The study of finite size scaling of a Taylor series presents a formidable technical task. Since 
the coefficients are generalised susceptibilities, each of them exhibits non-trivial finite size scaling. 
The scaling of the individual coefficients, evaluated at /I = 0, has to combine to the correct scal- 
ing of the finite density pressure given by the sum, thus requiring delicate cancellations between 
the individual contributions in the large volume limit. Classifications of the behaviour of various 



generalised susceptibilities are given in [34]. 



4.1 Quark number susceptibility to order /i for Nf = 2 

New results from this approach were reported this year by Gavai and Gupta $j2\\. They per- 
fomed simulations on L 3 x 4 lattices with L = 8 — 24, using the standard staggered action and the 
R-algorithm. The quark mass was fixed in physical units to m/To = 0.1. The aim of the simula- 
tions was to bracket the critical point by computing the Taylor coefficients of the quark number 
susceptibility up to sixth order (i.e. 8th order for the pressure) for various temperatures in the range 
T /To = 0.75 — 2.15, and extrapolate to finite fx. This was done for different lattice volumes in order 
to get an estimate of finite voulme effects. 

The results for the convergence radius Eq. (4.3) are shown in Fig. [7|. A rather strong vol- 



ume dependence is apparent. While for the smaller 8 3 lattice the estimators p n ,r n do not seem to 
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Figure 7: Estimators of the radius of convergence, Eq. (4.3), at T/Tq = 0.95. p„ (left) and r„ (right) vs. order 
n on 8 3 (red circles) and 24 3 (blue squares) [^2|. 



converge to a finite radius of convergence, the results on the larger 24 3 lattice are consistent with 
settling at a limiting value. The boundary between the two behaviours was observed to occur at 
Lm K « 5 — 6 or L « 16 — 18, with larger volumes tending to a smaller radius of convergence. It is 
reassuring that the numerical values for p„ and r„ are consistent with each other. Taking the large 
volume result at face value and extrapolating to all orders the estimate for the location of the critical 
point is [i$/T = 1. 1 ± 0.2 at T/T Q = 0.95, or \i%jT c = 1 .1 ± 0.2 [||]. 

4.2 The pressure to order \i 6 for Nf = 2 

Another investigation of the finite density phase diagram of the two-flavour theory was made 
by the Bielefeld-Swansea collaboration, also using the Taylor expansion of the pressure. This group 
works with a 16 3 x 4 lattice with p4-improved staggered fermions and a Symanzik-improved Wil- 
son action, simulating with the R-algorithm, the quark mass is set to m/7b ^ 0.4. The calculation 
to order /I 4 was performed in while new results on /I 6 are presented in [31 ]. The last work also 



contains detailed discussions of two interesting analytic calculations to compare with, namely the 
pressure in high temperature perturbation theory [j35|], which is going to hold at asymptotically high 
temperatures, as well as the hadron resonance gas model, which gives a rather good description of 
the pressure in the confined phase [|3"6|]. 



In agreement with [32] and qualitative expectations, their detailed results for the coefficients 
in the pressure series satisfy C(, <C C4 <C C2 for T > To, i.e. one would have coefficients of order 
one for an expansion in (fi/nT). An impression of the convergence of the series can be obtained 
by looking at the quark number susceptibility calculated to consecutive orders, as shown in Fig. ||. 
For T < 1.27b, the series seems to converge rapidly and the /x 6 -result is compatible with the one 
through order jit 4 . Around the transition temperature To, the /i 4 -results show a peak emerging with 
growing /x/Tq, which in [3C] was interpreted as evidence for a critical point. However, the /x 6 
contribution suggests that in this region results do not yet converge, and the structure is hence not 
a significant feature of the full pressure. 

Another analysis of the data in [ [3l| ] is devoted to a study of the convergence radius. Fig. |9] (left) 
shows the ratio of the Taylor coefficients at consecutive orders. In the relevant region T < To, the 
data do not seem to settle on a limit value. More strikingly, the data appear to fall right onto the solid 
lines marking the prediction for those ratios from the hadron resonance gas model, on both sides 
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Figure 8: Quark number susceptibility computed through (9(/i 4 ) (dashed lines) and (9(/i 6 ) (solid lines) [pip. 
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Figure 9: Left: Ratios of the expansion coefficients, Eq. (4.1 ). Right: Quark number density divided by the 
susceptibility. This observable goes to zero at a critical point. Solid lines in both plots represent predictions 
from the hadron resonance gas (T < Tq) and the Stefan-Boltzmann ideal gas (T > 7q). [|T]]. 



of the transition. The hadron resonance gas model does predict a Hagedorn-like deconfinement 
transition, but as a smooth crossover rather than a real phase transition. Thus the data do not give 
any indication of a critical point. This is corroborated by Fig. ^ (right), showing the quark number 
density normalised on the susceptibilty, n q /% q = J^-. This quantity is related to the compressibility 
in the plasma, and should go to zero at a second order phase transition point. 



Thus the conclusion in [ 3 1 ] is that there is no evidence for a citical point from these data. This 



conclusion is not in conflict with the results from [32] discussed in the previous section, since the 
simulations were done at a much larger quark mass, for which one would expect a critical point to 
be at larger values of fi. Moreover, a different action was used, making a direct comparison difficult. 
However, the conclusion is different from the earlier one by the same group based on \i results 
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[30], which were interpreted as showing evidence for a critical point. This highlights the need for a 
careful examination of as many terms in the series as possible, before results are conclusive. Indeed, 
previous experience with inferring phase structures from convergence properties of strong coupling 
series in spin models [|3|] shows that this is very much an "experimental science". Between 10- 
20 terms are known in some of these expansions, and while for some models stunningly good 
predictions about non-analytic behaviour are obtained, others still fail even at this high order. 

5. QCD at finite isospin density 

Another way to learn about QCD at finite density is by taking recourse to theories without a 
sign problem, which are sufficiently close enough to the physical situation of interest. I shall not go 
into the long list of activities along those lines, but concentrate on QCD at finite isospin density with 



chemical potential jtt/ [37]. For small enough Hi/T it can be argued that this theory should agree 



quantitatively with that at small jj,/T ', and recent lattice simulations support this picture [^, 39]. 

QCD at finite isospin density is obtained from two-flavour QCD by assigning opposite chemi- 
cal potentials to the quark flavours, \i u = —\l<i = fi, leading to \ii = (}X U — }Xd) = 2fi. This results in 
cancelling the phase of the determinant, so that the partition function now contains only its modulus 
and thus has real positive measure, which can be simulated without problems, 



DU\detM(n)pe-^ [U K (5.1) 

A schematic phase diagram of the theory is shown in Fig. [l(J On the lower right it features a pion 
superfluid phase, due to pion condensation ^ when |jU/| > m % . Note that in nature one 
cannot have a system with ^0 and ji B = 0, since the weak interactions do not conserve isospin. 
The interest in this theory is because of its formal relation to QCD at finite baryon density. It is 
also in this formal sense that the concept is generalised to Nf = 3. Indeed one would expect that 
for fa sufficiently small it should recover the physics at small baryon density, and hence the dashed 



transition line in Fig. 10 should be approximately the same as in the theory with baryon density. 
The argument goes as follows p8|]. An expectation value evaluated at finite baryon density can be 
rewritten as an expectation value evaluated at finite isospin density by means of the reweighting 
formula, 

Now consider probing the deconfinement transition with a gluonic observable O, e.g. the plaquette 
susceptibility showing a peak. As long as (cos 6}^ ~ 1, the observable O' = e' d O will signal the 
same transition as O. But in this regime one may as well neglect the phase altogether, in which case 
one probes the transition at finite isospin density. Based on this argument, one expects the transition 
lines in the two theories to be close to each other as long as reweighting works, i.e. for ju/7b < 1. 
(In an expansion about 6 = 0, the difference should be of the order ~ (0 2 ))- This expectation 
was numerically verified for the pseudo-critical surface 7o(m,/i) in the Nf = 2,3 theories. In 
the two-flavour theory, the Bielefeld-Swansea collaboration performed a Taylor expansion both 
in baryon and isospin chemical potential, and the resulting pseudo-critical lines were found to 
quantitatively agree [pOf]. Similarly, quantitative agreement was found between the pseudo-critical 
lines determined from finite [38, 39] and imaginary chemical potential /I,- [^, |6]]. 
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Figure 10: Schematic phase diagram for QCD at finite isosipin density. 



6. Systematics of reweighting and Taylor expansion 

A few years into their existence, there are now several investigations of the systematics of the 
reweighting and Taylor expansion approaches. Recalling the double reweighting formula 

(0) iP a) = ^ >(M) ~ <T VAF , (6.1) 

one would like to estimate when the exponential suppression of the signal becomes insurmount- 
able. Splitting the determinant into modulus and phase, detM = | detM|e' e , this should occur when 
(cos 8) <C 1, or equivalently when the root of the variance of the phase of the determinant grows 
larger than n/2, 



a(6) = ^<0 2 ) - (d) 2 = ^(d 2 ) > n/2. (6.2) 

In order to quantify this, the Bielefeld-Swansea collaboration evaluated the phase by means of its 
Taylor expansion [J3l|], 

.00 = *___£ ^ (63 ) 



Contours of values for the variance are shown in Fig. |11|. According to the criterion Eq. (p.2|), the 
line corresponding to a = n/2 can be viewed as the boundary for the reliability of reweighting. The 
region to its lower right is safe while one would not trust results obtained up and left from it. This 
means the deconfined phase of QCD is rather accessible as expected, while in the transition region 
one finds the constraint {J. /To < 1, in accord with the constraint on other methods. The figure shows 
contours for one given volume. Reweighting gets exponentially harder with volume, and thus the 
contour lines move rapidly to the lower right as the volume is increased. 



In another paper [41] Ejiri discusses the difficulties of a combined application of reweighting 



and an analysis of Lee- Yang zeros (LYZ) [21]. The latter exploits the fact that on a finite volume 
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T/T„ 



Figure 11: Contours of a (9) for a fixed volume 16 3 [pT|], 



there are no singularities in the pressure, and hence no zeroes of the partition function for real 
couplings j8 . However, there are zeroes for comlex couplings, whose real part indicates the location 
of an analytic crossover, i.e. the pseudo-critical temperature Tq. In the infinite volume limit, these 
zeroes move to the real axis if there are true phase transitions, while they stay at complex values 
for crossovers. A LYZ analysis for reweighted finite /I then amounts to numerically searching for 
zeros in the expression 



2 (j3 Re ,j3i m ,A0 



,(JVf/4)(lndetM(/t)-lndetilf(0)) 



(6.4) 



where one additonally reweights into the complex coupling plane. Ejiri argues in [[fj] that this 
combined procedure does not have an infinite volume limit. Taking V — > °° at finite statistics, the 
above expression for the partition function will always go to zero because of the sign problem, and 
hence always signal a phase transition, even where there is a crossover. This point of principle 
is not surprising. Indeed the same mechanism precludes a numerical infinite volume limit of any 
observable computed via reweighting. However, the question for practical simulations is whether 
for a given volume enough statistics can be gathered to beat the sign problem, and whether the 
volume is large enough to reproduce infinite volume physics with sufficient accuracy. Eq. Q5.4| ) 
illustrates the difficulty of this procedure: the LYZ get masked by the noise from the reweighting 
factor, and one has to guard against mistaking a disappearing signal for a Lee- Yang zero. The 
problem boils down to being able to give reliable errors for the reweighting procedure, which are 
needed for a qualified judgement on whether statistics is sufficient or not. 

In an interesting qualitative investigation of systematics, Splittorff makes use of the finite 



density formulation [|42|]. He suggests to turn the reweighting argument for approximate equality 
of finite isospin and baryon density around, in order to determine the limit of applicability for 
reweighting. For this purpose a matrix model prediction for the transition line to the pion 
liquid is combined with the contour lines for the variance of the phase of the determinant, Fig. 1 1 



as shown in Fig 12 (left). The value of a (6) rises towards the lower right, and one observes that 
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Figure 12: Left: rescaled version of the contours of <j(9), Fig. [flj The values of the contours increase 
towards the lower right and approach the transition line to the pion condensate at finite isospin density. 
Right: the critical endpoints for two sets of quark masses by the reweighting method [^0| ^J. From [|42|]. 



contours become denser with larger values, approaching the transition line to the pion liquid. This 
is to be expected on the grounds that a non-trivial phase of the determinant wipes out the pion 
condensate, which is not present at finite baryon density. Hence Splittorff suggests to interpret 
the transition line to the pion liquid in the theory at finite isospin density as a "cut-off" for the 
applicability of reweighting: continued reweighting to the right of that separating line would mean 
one has a serious overlap problem, since the presence of the phase makes a physical difference 
there. That this is a matter worth exploring in more detail is shown in Fig. |l2| (right), which 
displays the critical endpoints from reweighting for two different quark mass sets [20, g], and both 
fall in the neighbourhood of this boundary. 

A similar argument may be applied to results from the Taylor expansion. For larger volumes, 
the sign problem becomes more severe. In the Taylor expansion, whose coefficients are evaluated 
at jJ. = 0, this shows up in two ways. Firstly the need for more terms to describe the sharpening 
divergence in its build-up. Secondly, the need for ever more precise cancellations between different 
terms in order to combine to the correct volume scaling behaviour of the sum. But the severity of 
the sign problem also puts a limit on /I /T for fixed volume, as we have seen already. Checking in 
Fig. [8] at which values of /I /7b the sixth order contribution to the susceptibility becomes important 
for a given temperature, Splittorff concludes that the 4th order expansion only works to the left of 
the leftmost contour line in Fig. 12 (left). Approaches based on imaginary /x never face the sign 
problem. However, as dicussed in the next section, analytic continuation to real /x necessitates a 
Taylor expansion too, and one would expect a similar limitation. Of course, these estimates are not 
yet quantitative, as the finite isospin transition line is determined from a model and not known with 
any accuracy, but they point out interesting directions to pursue. 



7. QCD at imaginary ji 



Since the QCD fermion determinant with imaginary fx = i/Xj is real positive, it can be simulated 
just as for jx = 0. It is then natural to ask whether such simulations can be exploited to learn 
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H/(tcT) H B /T 

Figure 13: Left: Schematic phase diagram for QCD at imaginary chemical potential. Right: Nf = 3 results 
after continuation |Q and compared to Nf = 2 + 1 [p0|]. The light quark masses are the same. 



something about physics at real }X. The strategy to get back to real jj, is to fit the Monte Carlo 
results, wich are free of approximations, to a Taylor series mjx/T. In case of apparent convergence 
it is then easy to analytically continue the power series to real }X. This idea was first used for 
observables like the chiral condensate and screening masses in the deconfined phase [ J44| , [27| ]. It was 



then shown to be applicable to the phase transition itself [|40|], which has recently been exploited in 
a growing number of works ||6[]-p9|]. 
The partition function, 

Z{V,pLj)=TvU-^-^l T ), (7.1) 



is periodic in the imaginary direction, and the period can be shown to be 2n/N c for N c colours 
[p0[]. Hence, in addition to being even in pL, the QCD partition function has the additional exact 
symmetry Z(fi r /T,fii/T) = Z(fX r /T \pii/T + 2n /3). Because of the fermionic boundary conditions 
in the Euclidean time direction, this symmetry implies that a shift in /I, by certain critical values 
is equivalent to a transformation by the Z(3) centre of the gauge group. Thus, there are Z(3) 
transitions between neighbouring centre sectors for all {}Xi/T) c = ^ (n+ j) ,n = 0,±1,±2,.... 
It has been numerically verified that these transitions are first order for high temperatures and a 



smooth crossover for low temperatures [40, 45]. As a consequence, the schematic (T,Hi) phase 



diagram looks as in Fig. |13|. The vertical line coming from the top denotes the Z(3) transition, 
while the deconfinement transition line now bends upwards as a function of /i,-. The order of the 
transition and the existence of an endpoint depends again on the number of flavours and the quark 
masses. Because of the symmetry of the partition function this picture is then periodically repeated 
for larger values of Hi. 

The idea then is to simulate with imaginary chemical potential, and fit the full simulation 
results by a power series of order N, 

N /li-\2n 
n 

Since the Monte Carlo results contain no approximation or truncation, convergence can be in- 
spected by the quality of the fits to the data. In the case of satisfactory convergence analytic 
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continuation fa — > ifr is a trivial matter. It was shown that this strategy can be extended to 
the pseudo-critical line itself, which on finite volumes is a smooth function with an even Taylor 
expansion [Boh. Detailed comparisons give quantitative agreement for To(m,/i) computed from 



imaginary jj. and other methods []6} |45|], cf. Fig. |13[ There are also first results for Nf = 4 with 



Wilson fermions [49] 



This procedure may be expected to converge as long as the value of /x, does not exceed the 
critical value of the first Z(3) transition, \fl\/T < It/3, or /Xe^SSOMeV in physical units. This 
constraint is due to two limitations. Firstly, in the infinite volume limit the location of the Z(3) 
transitions would bound the radius of convergence of the Taylor series. And secondly even on 
finite volumes, where there are no non-analyticities, no new information is obtained by going to 
larger /I,- because of the periodicity. Nevertheless, interesting arguments are being made that one 
might well extend the continued results along the real /i-axis beyond this radius with the help of, 



e.g., Pade approximants jplQ. 

Working at imaginary /I has a couple of technical advantages. It is computationally simple and 
much cheaper than reweighting or computing coefficients of the Taylor expansion. Moreover, both 
parameters j3 , jJ. are varied and thus one obtains information from statistically independent ensem- 
bles. It also offers some control on the systematics by allowing a judgement on the convergence 
of the fits. Furthermore, it is a good testing ground for effective QCD models: analytic results can 
always be continued to imaginary pL and be compared with the numerics there, as demonstrated for 
several examples in [H^J. The main limitation presently is the radius of convergence in the large 
volume limit, /l/T ~ 1. 

7.1 A generalised imaginary ji. approach 

The method of simulating imaginary /I and analytically continuing can be generalised in an 



interesting way, as suggested by Azcoiti et al. [|47|]. The idea is to rewrite the standard expression 



for the staggered fermion action at finite density Q52[ ] by replacing the chemical potential with two 
new parameters x,y, 

1 



~* x \ £ YnVoin) (U n fiYn+0 ~ U l-Q.oYn-o) +y\ £ 9nHo(n) (u n flYn+0 + V«-o) (7.3) 

where x = cosh(a/x), y = sinh(a/x). This means the action has been enlarged by an extra parameter. 
The ordinary finite density action is recovered by the constraint x 2 — y 2 = 1. Thus, if the solid line in 



Fig. [14j denotes a phase transition line in the x,y plane of the enlarged theory, its intersections with 
the dotted line representing the constraint correspond to physical transition points. The enlarged 
theory still has the sign problem, but one can simulate at imaginary y = iy. The potential of this 
method to improve over the simple imaginary jit approach is that there are now different parameter 
sets x,y to be simulated, so one might hope to be able to extrapolate in a controlled way to reach 
larger values of pL /T and thus probe the phase diagram at lower temperatures. 

Numerical results from this method for the four-flavour theory have been presented in 
Simulations were performed on 8 3 x 4 lattices with standard staggered fermions and the R-algorithm. 



Fig. 15 shows results for the pseudo-critical coupling obtained at imaginary y, fitted to the form 
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Figure 14: Schematic phase diagram for QCD at imaginary chemical potential. 
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Figure 15: The critical coupling in the Nf — 4 theory as a function of imaginary y. The leading order 
quadratic fits are continued to real y, the intersection with the vertical line gives the physical value [|48|]. 



Pc(y) = fk> + jSiJ 2 and continued to real y. The vertical dotted line gives the intersection with phys- 
ical QCD. Clearly, the method works as well as the ordinary imaginary /I approach. However, its 



full potential of simulating different sets of x,y for a fixed }i/T has not been probed yet. Fig. |15 
(right) shows an example of a result that was extrapolated beyond the first Z(3)-transition, which 
shows up by the kink in the data. This is done on the grounds that no such kink is present in the 
real direction. Nevertheless, even for imaginary y, beyond the kink the curve is not constrained by 
any additional data points, and hence convergence of the extrapolation is not guaranteed. 

New and as yet unpublished results for the two-flavour theory with indications of a critical 



point are presented in these proceedings [53]. An interesting open question is whether the method 
can indeed be used to obtain more control over analytic continuation. Another promising idea by 
the authors is to use their action for reweighting. The two parameters x, y might be tuned such as 
to shorten the reweighting distance to the physical point of interest. 

7.2 Imaginary pL and Fourier transformation: QCD at fixed baryon number 

Last year promising attempts of an alternative use of imaginary /I have been made [|4|, |55| ]. 
This approach makes use of the relation between the grand canonical partition function at imag- 
inary chemical potential and the canonical partition function at fixed baryon number via Fourier 
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transformation pO, ES\. After earlier attempts on a Hubbard model [57], there are now promising 



new results on QCD presented in these proceedings [|58|, |59j]. The difficulty in this case is to make 
baryon number large enough so as to reproduce the finite chemical potential calculations in the 
thermodynamic limit. 

Baryon number is fixed by inserting 8(3B — J d 3 x YYoV) i nto tne P at h integral. The result 
is the canonical partition function, which is related to the grand canonical partition function at 
imaginary pL via a Fourier transform, 

Z C (B) = -L f d (£) e- BB TZ GC (ii = fa). (7.4) 



2% J -n V T 

The idea followed in []54j] is to sample Zgc(A 4 = i^Mc) by Monte Carlo, and then compute 
Z C (B) I 1 



■ J dUi exp (i3B^) det(irt) \ , (7.5) 



Zac{iHMc) \det(//i MC ) 

i.e. the fermion determinant gets Fourier transformed, and so it has to be calculated exactly. This 
is costly, but the benefit is that now no analytic continuation or Taylor expansion is needed. The 
sign problem of course resurfaces here as well, making Zc(B) noisy. But the strength of this effect 
is governed by baryon number, and not by volume directly. In the thermodynamic limit one would 
have to send baryon number to infinity in order to have fixed baryon number density. Thus, the 
larger the volume, the smaller the accessible baryon number density for a simulation. Nevertheless, 
fixing a small baryon number makes sense in order to study e.g. nuclear few body systems, and as 
the following results show, with sufficient computer power one can reach reasonably large baryon 
numbers to make contact with the grand canonical formulation. 

Numerical results obtained by de Forcrand and Kratochvila are shown in Fig. [l^. They were 
obtained on a 6 3 x 4 lattice with Nf = 4 standard staggered fermions and hybrid Monte Carlo 
simulations, the quark mass was m/To sa 0.2. The left panel shows the conversion from baryon 
number to chemical potential. This is achieved by evaluating the free energy F(B) as function of 
different fixed baryon numbers, and computing the grand canonical partition function by Laplace 
transformation. The integral can be evaluated by means of a saddle point expansion, 

Z GC (At) = Jdp^(~(f{p)-fip)\ (7.6) 

yielding the chemical potential as a function of baryon number, 

,, , F(B + 1)-F(B) 
M~/(P)~- { — ■ (7-7) 



As Fig. 16 (left) shows, on a 6 lattice it is possible to perform this up to quite respectable 
baryon numbers. The resulting picture essentially shows the Maxwell construction for changing 
between canonical and grand canonical ensembles. The S -shaped curves are indicative of a first 
order phase transition and represent the metastability region. The two envelopes can be ascribed to 
the hadronic and quark gluon phases and are well fitted by hadron and weakly interacting massless 
gas models. Preliminary results from the canonical ensemble obtained with Wilson fermions have 
been published in [|55|] and are presented in these proceedings 
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Figure 16: Left: Numerical Maxwell construction to relate the canonical and grand canonical ensembles. 
Right: Comparison of different methods. Agreement is quantitative for jx/T<l. From [p8H . 



Fig. 16 (right) shows the resulting critical line for the four flavour theory in comparison with 
ones obtained by other methods discussed here. All data are generated with the same action and 
for the same quark mass and lattice spacing, only the volumes differ between the data sets. For 
H/T<1, the quantitative agreement is impressive. Only the data from [ fl6| ] are somewhat off the 
others, presumably due to the much larger volume of that data set. Note that for n/T > 1.3 agree- 
ment stops and the different data sets diverge. This is in accord with the previous statements that all 
methods discussed have roughly the same range of applicability, but different systematics. The data 
continued from imaginary /I in this region are essentially unconstrained and just extrapolated. For 
the reweighted data points the expectation value of cos 6 is quoted in the figure for selected points. 
For /J. /T > 1.2 it is completely lost in the noise and hence the data points are not trustworthy. 
Note also that the data coming from the canonical ensemble in principle do not have the restriction 
/a/T<1. The data bend down more strongly as one might expect in this region of the phase dia- 
gram. It will be exciting to see whether the curve can be reliably continued beyond ~ 1 . With 
the density of states method [24], Pade approximants [51] and the canonical ensemble [Q|, there 
are at least three attempts at work in this direction. 



8. The critical end point and its quark mass dependence for three flavours 

As has become clear by now, a determination of the order of the phase transition and the critical 
point is much more demanding than the location of the pseudo-critical temperature 7b(ra,//). The 
best starting point for such an enterprise is the Nf = 3 theory. This is because we know there is 
a critical point at ji = 0. In the phase diagram Fig. || (left), the critical quark mass m c {ji. = 0) 
separating the crossover from the first order sections along the three flavour diagonal is known to 
be at a moderately small value accessible to simulations [Q, g, £0. With the quark mass tuned to this 
value, the first order transition line in the (T,/j. 2 ) phase diagram Fig. 17 (left) reaches all the way to 
the temperature axis, on which it ends with a critical endpoint. According to the standard scenario 
discussed in the introduction, if we now increase the quark mass to values m c > m c (0), the whole 
transition line shifts, with the critical endpoint wandering to the right towards a real jx c ^ 0, thus 
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Figure 17: Left: Schematic phase diagram in the (T, /z 2 )-plane for different quark masses. Solid lines are 
first order, dotted lines crossover, and the thick line represents the endpoints. Right: The critical quark mass 
as a function of fi 2 . The vertical line on the left marks the first Z(3) transition at imatinary fi. From [ph. 



tracing out a smooth function fl c (m). On the other hand, if we change the quark mass to m < m c (0), 



the critical endpoint should wander in the imaginary \x -direction to the left, as shown in Fig. [17. 

Inverting the function fl c (m), we are interested how the critical quark mass separating first 
order from crossover changes with pL. This function again has a Taylor expansion in even powers 
of }x, and one expects 



+ ... 



(8.1) 



m c (n =0) 

with coefficients of order one. In the three-dimensional phase diagram of Fig. [3] (right) this means 
we are looking for the curvature of the critical surface in the Nf = 3 direction at w c (0). Once this 
functional dependence is determined, it will return the critical end point pL c for a given quark mass. 

8.1 Numerical results for Nf = 3 from imaginary }X 

This program was recently carried out in using 8 3 x 4 lattices with the standard staggered 
action and the R-algorifhm. The observable used to find the critical point was the Binder cumulant, 
making use of the knowledge that the transition is in the universality class of the 3d Ising model. 
For a critical point in this universality class the value of B4 is accurately known, 



B 4 (m c ,n c 



((SW) 4 ) 
((SW) 2 ) 2 



1.604, V^oo, 



(8.2) 



while B4 — > 1(3) for a first order transition (crossover). Hence in the infinite volume limit B4 
is a non-analytic step function. However, on finite volumes it will pass through the Ising value 
smoothly, with a slope increasing with volume. Measurements for several values of imaginary 



/i and quark masses are shown in Fig. 18 (left). The data can be fitted to a leading order Taylor 
expansion in both the quark mass and chemical potential about the known critical point at m c {\i = 
0), 



B 4 (am,aiJ.) = 1.604 + B (am — am c (0) + A(aju) ) + . . 



(8.3) 



From the fit parameters one can directly extract the desired coefficient d(am) / d(aji) 2 in lattice 
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Figure 18: Left: The Binder cumulant calculated for various masses and chemical potentials. The value 
1.604 corresponds to a critical point. Right: Finite volume scaling ~ L 1//y of the fit parameter B, Eq. (S.3). 
v = 0.62(3) is consistent with the universal Ising exponent v = 0.63 [0]. 



units. For the continuum conversion one needs to take into account that the lattice spacing effec- 



tively changes with r(ju), and hence a(ju) ^ a(0). For c\ in Eq. (8J_) this yields 



1 dm c 7l 2 d(am c ) 1 dT 

Cl = m f (0) d{\Ll%Tf = N?(am c )(0) d{a^i) 2 + %d{ii/nTf (-) 



Fig. [18| (right) shows a finite size scaling analysis of the fit coefficient B, and the fitted volume 
scaling nicely reproduces the exponent predicted by universality even on moderate volumes. Note 
that these are extremely difficult calculations, as very long Monte Carlo trajectories of order 80k 
are required in order to get sufficient tunneling statistics. 

In accord with qualitative expectations, c\ m 0.8(4) ^ is of order one. However, the large 
statistical error indicates that the coefficient is also consistent with being close to zero. Independent 
of the accuracy of the result, an important conclusion is that m c changes very little when /I is 
switched on, or conversely that pL c changes rapidly under small variations of the quark mass. 

8.2 Numerical results for Nf = 3 at finite isosipin density 

In a recent article Kogut and Sinclair report on similar investigations in the theory at finite 
isospin density [p9[]. They work on L 3 x 4 lattices with L = 8 — 16, using the standard staggered 



action and the R-algorithm. Their quark mass is chosen as m>m c (0), with the aim to see how the 
critical point moves as a function of /I/. As an observable they use the Binder cumulant discussed 



in the previous section. Results from their simulations are shown in Fig. |19[ The left panel displays 
an investigation of step-size effects on the Binder cumulant, which are found to be significant in this 
quark mass regime. Instead of extrapolating to zero step size, the standard usage of the R-algorithm 
is to simulate at some reference step size whose error is known to be smaller than typical statistical 
errors at some reference mass m. When going to smaller quark masses, a common practice is to 
keep the step size to be half the bare quark mass. However, in the low quark mass regime of interest 
here, this procedure clearly breaks down, with even qualitative changes of the results. The figure 
shows how at the same quark mass the transition looks clearly first order for large step sizes, but 
changes to crossover behaviour once extrapolated to zero stepsize, which thus is mandatory. 
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Figure 19: Left: Stepsize dependence of the Binder cumulant. Right: Binder cumulant extrapolated to zero 



stepsize. Dependence on isospin is very weak [ 39 1 



The right panel shows the results after such extrapolations have been performed. In agreement 
with the findings from imaginary /I in the previous section, the Binder cumulant is practically flat 
and only very weakly depending on Hj. All data points are well in the crossover regime. What is 
surprising is that the weak /^-dependence has a tendency for dB^/dpLj > 0. This would imply that 
as fa is switched on the transition moves deeper into the crossover regime, instead of approaching 
a critical point! 

8.3 Nf = 3 and Nf = 2 + 1 with an exact algorithm 

In order to check for stepsize effects and clarify the sign of the /I -dependence of B4, de For- 
crand and I have redone the Nf = 3 calculation of ^ reported in Section [O] with the exact rational 
hybrid Monte Carlo (RHMC) algorithm developed by Clark, Kennedy and Sroczynski [|60|] (also 



presented in these proceedings [61]). In a first test we compared simulations of the Binder cumu- 



lants performed with that algorithm to ones from the R-algorithm extrapolated to zero stepsize. The 



results are shown in Fig. 20 (left), and perfect agreement is found. The right panel then shows a 
new determination of the critical quark mass m c (ju), both for zero density and an imaginary chemi- 
cal potential, using the RHMC algorithm. Note that the Binder cumulant now passes its Ising value 
at a significantly different mass compared to the results in the literature. We find am(0) sa 0.026, 
which is a shift of about 25% due to stepsize effects! On the other hand, switching on a chemical 
potential has no effect on the Binder cumulant. As a preliminary result we find d{am c ) / d{a}i) 2 « 
within errors, which is consistent with the findings at finite isospin [[3^] reported in the last section. 



We have also mapped out the critical line for non-degenerate quark masses, as shown in Fig |2 1 
both with the R-algorithm and the RHMC algorithm. In general there is a significant step size effect. 
The picture clearly puts the physical point, where Fodor and Katz performed their simulations, on 
the crossover side of the line. Note, that the physical point is very close to the critical line. This 
is consistent with the requirement of finely tuned quark masses in order to have a critical point at 
moderate chemical potentials. (The calculation of c\ in this case is still in progress. Taking our 
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Figure 20: Left: Comparison of the Binder cumulant computed with the RHMC algorithm (leftmost data) 
and the zero stepsize extrapolation of the R-algorithm. Right: Determination of m c (jj,) for ji = 0,0.2 with 
the RHMC algorithm. The arrow marks the result from the R-algorithm. 
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Figure 21: The critical line separating first order from crossover for Nf = 2 + 1. A significant shift is 
observed when eliminating stepsize effects. The line passes in close vicinity of the physical point (FK). 



R-algorithm result c\ from the three flavour case, our resulting \i c would be consistent with theirs 
within errors). If the chiral limit of the two flavour theory turns out to be 0(4), there is a tri-critical 
point at some quark mass m[ nc on the m r axis. Our results can be fitted with the corresponding 
scaling equation and would then predict nf™ /Tq m 2.8. 

8.4 A non-standard scenario for the phase diagram 

Let us assess the consequences of the step size effects on the critical point. After continuum 
conversion the new result for m c ( i u), now free of step size errors, is 



1-0.6(2) f-^) 2 . (8.5) 



m c (n) 

m c .(ju=0) " "'" V_/ V7rry 

Note that the sign of the leading term has changed compared to the previous result! The reason is 
that the first term in Eq. ( |8~4] ) now is consistent with zero, so the negative second term dominates. 
This means the critical mass gets smaller when a real ji is switched on, and hence that the critical 



surface in the phase diagram leans towards smaller quark masses, Fig. 22, i.e. the opposite of the 
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Figure 22: For dm c (jJ.)/dji 2 < 0, there is no critical point at all, the dotted line on the right is merely a 
crossover. 



standard scenario Fig. [3] (right). In this case, the first order region in a plane of constant /J. is actually 
shrinking for growing /x. If the physical point is in the crossover region at /I = 0, then switching 
on a chemical potential will not lead to an intersection with the critical surface, and hence there 
would be no critical point or first order phase transition at all! Note that this scenario is perfectly 
consistent with all the universality arguments summarised in Section |Q| . The (r,ju) phase diagram 
would then only have the transition line separating the superconducting phase from nuclear matter, 



as in Fig. 22 



8.5 Can one expect a critical end point at Hb ^ 500 MeV? 

Clearly, as the discussion of the systematics in the last section revealed, this last result is pre- 
liminary as well, an important question being in which direction corrections go as the continuum 
limit is approached. Nevertheless, based on the results on the curvature of the critical surface, 
one may obtain a rough estimate for the conditions required to have a critical endpoint in the phe- 
nomenologically interesting region /j.b < 500 MeV. Irrespective of whether an eventual continuum 
result for m c (jj.) will have positive or negative curvature, as long as the coefficient ci in Eq. ( |rO| ) is 



~ 0(1) ( its natural size; all known coefficients in the pressure [31], screening masses p7[ ] and the 
pseudo-critical temperature [^] are of that order), it implies a very strong quark mass dependence 
of the value of /i c . For instance, in order to have /i c ~ 120 MeV as predicted by Fodor and Katz [||] 
for Nf = 2 + 1 , the quark mass has to obey the condition 

K— t^-^<1.05, (8.6) 
m c (n =0) 

i.e. it has to be fine-tuned to be within 5% of the critical quark mass. Provided the coefficient c\ 
does not change drastically in the case of Nf = 2 + 1 , a similar situation will be encountered there 
as well. 

At this point it should become clear that we are still far from a quantitative solution of the 
problem of the critical endpoint. To achieve a resolution better than 5% in the quark masses would 
even require to distinguish the up and down quarks. By contrast, we have just discovered a 25% 
systematic error in the critical quark mass due to step size effects, and we have discussed earlier 
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Figure 23: The critical heavy quark mass separating first order from crossover as a function of pi 1 



the strong cut-off dependence (~ 100%) of the critical quark masses in physical units. Thus we 
should expect formidable shifts in jx c on the way to a reliable continuum result. 

8.6 The heavy quark limit: Potts model 

Other interesting projects are concerned with the upper right corner of Fig. [3], i.e. the region 
towards the quenched limit. Simulations of quenched QCD at finite baryon number have been 



done in [63]. As the quark mass goes to infinity, quarks can be integrated out and QCD reduces to 
a gauge theory of Polyakov lines. First simulations of this theory with Wilson valence quarks can 
be found in [|64j]. At a second order phase transition, universality allows us to neglect the details 
of gauge degrees of freedom, so the theory should be in the universality class of the 3d three-state 
Potts model, which is the 3d Ising model. Hence, studying the three-state Potts model should teach 
us about the behaviour of QCD in the neighbourhood of the critical line separating the quenched 
first order region from the crossover region. For large pL the sign problem in this theory was actually 
solved by means of cluster algorithms recently [|6^]. 

Here I want to discuss an as yet unpublished result on simulations of the three state Potts 
model as presented in these proceedings []65[|. For small fl/T, the sign problem of this theory is 
mild enough so that brute force simulations at real pL are feasible. In the simulations presented 
in j65|], the change of the critical heavy quark mass is determined as a function of real as well 



as imaginary \x, as shown in Fig. [23[ Note that M c {}i) rises with real chemical potential, i.e. the 
first order region in Fig. |5| shrinks as finite baryon density is switched on. This system is thus 
an example of the non-standard scenario discussed in the previous sections! Note also that the 
qualitative behaviour in going from real to imaginary }x is exactly as predicted in the schematic 
picture Fig. [l7| [[J, and analytical continuation in determining this critical line thus works. 



9. Conclusions 



The last couple of years have seen an enormous increase in activities concerned with lattice 
determinations of the QCD phase diagram in all of its interesting regions and limits. While definite 
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conclusions cannot yet been drawn, there is a lot of progress in refining the methods and studying 
the systematics. 

The longstanding question of the nature of the phase transition in the two flavour theory in the 
chiral limit is still open. But large volumes are now available, and simulations on N t = 6 lattices 
will be undertaken soon. In combination with now available exact algorithms these will hopefully 
settle the issue in the near future. 

There are now several groups that are tackling the critical endpoint. However, these investi- 
gations are extremely difficult and still carried out over a scatter of theories and parameter values. 
Within the R-algorithm, the critical endpoint from two-parameter reweighting is consistent with 
the shape of the critical surface determined from imaginary chemical potential. However, the R- 
algorithm in the regime of physical quark mass has been demonstrated to be afflicted by strong 
stepsize effects, which change the apparent order of the phase transition. Exact algorithms are now 
being employed successfully, and this source of error will soon be eliminated. An important qual- 
itative conclusion is that the critical chemical potential of the endpoint is extremely quark mass 
sensitive. A critical point <400 MeV requires the physical light quark masses to be less than 
5% larger than the critical values at zero density. While it is quite possible that nature has arranged 
for this, it is clear that under those circumstances a quantitative determination is going to be a 
formidable task: any systemtic error in the current simulations is going to have enormous effects 
on the location of the critical point. Recall that all calculations reported here are on coarse lattices 
with a ~ 0.3 fm, and in most works quark masses are only fixed in lattice units. Furthermore, 
finite volume and stepsize effects have been shown to be larger than several 10%. Under those 
circumstances it is still conceivable that there is no critical point and phase transition at all. This 
means working towards producing results in the thermodynamic and continuum limits will be just 
as exciting as the first qualitative calculations! 

Acknowledgements: I am grateful to Philippe de Forcrand for valuable discussions and comments, 
help with figures, and a continued enjoyable collaboration. 
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